Configurational entropy of network-forming materials 
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We present a computationally efficient method to calculate the configurational entropy of network- 
forming materials. The method requires only the atomic coordinates and bonds of a single well- 
relaxed configuration. This is in contrast to the multiple simulations that are required for other 
methods to determine entropy, such as thermodynamic integration. We use our method to obtain 
the configurational entropy of well-relaxed networks of amorphous silicon and vitreous silica. For 
these materials we find configurational entropies of 1.02fcs and 0.97fc_g per silicon atom, respectively. 

PACS numbers: 71.55. Jv,61. 43. Dq,61.43.Bn 
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In materials such as vitreous silica, amorphous silicon 
or vitreous ice, the structure is determined by the set of 
bonds (covalent or hydrogen) between particles. In this 
manuscript, we will refer to these materials as network- 
forming materials. While the local environment of each 
particle is usually governed by strict rules, the bonded 
network can show a wide variety of different topologies. 
The focus of this work is to present a method to estimate 
the number of topologies that a network-forming material 
can take, or more precisely, its configurational entropy. 

The common computational procedure to estimate the 
entropy S at temperature T2 is to measure the energy E 
as a function of temperature and then to integrate from 
a temperature Ti at which the entropy is known: 



S{T2) - S{T,) - 



(1) 



This requires sampling a large number of different net- 
works and can therefore only be applied to systems with 
fast dynamics. In this work, we present a completely 
different approach to estimate the entropy, based on in- 
formation theoryld, and related to the work of Schlijpcr 
et al. who determiujad the entropy of the Ising and three- 
states Potts modelsa. 

One important concept in information theory is the 
Shannon entropy. It is commonly explained in the con- 
text of a string of n bits. In this case, the Shannon 
entropy H{n) is defined as: 



H{n) = -^p(i)log2p(i 



(2) 



where the index i runs over all possible bit sequences of 
length n and p{i) is the probability of sequence i occur- 
ring. A related concept is the entropy density s of a large 
string of N bits: 



s = lim [H{n+1) - H{n)] . 



(3) 



A practical procedure to estimate s for such a long string 
is to extract from it a large number m sub-sequences. 



each containing n bits with n <^ N. An estimate for the 
probabilities p{i) is then given by: 



p(i) « — , 

TO 



(4) 



where fi equals the number of times sub-sequence i was 
observed. The estimates for p{i) in combination with 
Eq. yield H{n}. The entropy density s is then ob- 
tained using Eq. {m . Usually, the limit n — > cxd converges 
rapidly and even moderate values of n are sufficient to 
predict s accurately. 

For systems in equilibrium it is easily shown that 
the Shannon entropy and the thermodynamic entropy 
of Eq. (|^) are equivalent, apart from a factor of ln(2). 
In this case, the probabilities p{i) are simply the Boltz- 
mann weights: p{i) = exp(— /3£'i)/Z, where /? is the in- 
verse temperature, Z the partition function and Ei the 
energy of state i. In the present work we show how infor- 
mation theory can also be used to obtain the configura- 
tional entropy of network forming materials. The method 
we present can be applied to any network provided (1) 
the atomic coordinates are known and (2) a list of bonds 
is supplied or can be constructed — for instance based 
on a distance criterion — which uniquely determines the 
network. 

To determine the configurational entropy of a network- 
forming material we choose a large number to of random 
positions in the simulation cell. For each position we find 
the nearest n particles and identify the graph formed by 
the bonds connecting these particles. We then assign_a 
label to this graph, based on the graph automorphismEI. 
We count the number of times a graph is observed and 
use Eq. (^) to estimate its corresponding probability of 
occurrence. These probabilities are fed into Eq. (|^) to 
obtain H{n). 

A small displacement of one of these random positions 
will usually result in exactly the same list of n nearest 
particles and thus the same graph. Consequently, there 
is an upper bound to the number to of random posi- 
tions that one should choose in a simulation cell con- 
taining N particles. An estimate of this upper bound 
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FIG. 1: Typical behavior of the corrected Shannon entropy 
Hc{n) in units of ks as a function of graphsize n; crosses (cir- 
cles) mark measurements for which H\{n) contributes less 
(more) than one percent to H{n). The dashed line is a 
straight-line fit to the crosses, starting from n = 10. The slope 
of this line is our estimate for the configurational entropy per 
atom a. These data are obtained from a two-dimensional sil- 
lium configuration containing 20,000 atoms. 



is obtained from the typical distance over which a ran- 
dom position can be displaced without altering the se- 
lected graph. From this we obtained as upper bounds 
m = 1.6nN and m = SAir^N for two-dimensional and 
three-dimensional networks, respectively. 

A related side-effect of choosing random positions is 
that the number of different graphs observed gets multi- 
plied with a factor proportional to n'^^^. This results in a 
correction to H{n) of the form g{n) = {d— 1) log(n), with 
d the spatial dimension of the network. The latter can be 
verified in crystalline networks where the configurational 
entropy is zero. 

Graphs with a probability p{i) smaller than 1/m will 
likely be observed only once, if at all. This finite-size 
effect grows with n, when the selected graphs become 
very complex. To monitor the impact of this effect, 
we record the quantity Hi{n), defined as the contribu- 
tion to H{n) of the topologies observed once. We reject 
the measurements for which Hi(n) exceeds one percent 
of H{n). Furthermore, we observed that the quantity 
H{n) + c ■ Hi{n) / H{n) converges much faster with in- 
creasing m than H{n) itself, with a suitable choice of 
c. We therefore use this extrapolated value as our best 
estimate for H{n). 

In Fig. |l| we show the typical behavior of the corrected 
Shannon entropy Hc{n) = H{n) — g{n) as a function of 
graphsize n. The entropy follows from Eq. (||) . The rapid 
convergence of the limit is demonstrated by the linear 
behavior of Hc{n) for intermediate n shown in Fig. |l| 
by the dashed line. The entropy per atom equals the 
slope of this line. In Fig. |^ we have marked with circles 
the points for which Hi{n) exceeds one percent of H(n). 
These points suffer from finite-size effects and should not 
be used. 




FIG. 2: A three-dimensional sillium network. Each particle 
is four-fold coordinated but no long-range order exists. 




FIG. 3: A two-dimensional sillium network. Each particle 
is three-fold coordinated, with a preferred bond angle of 120 
degrees. 



To test the usefulness ofj-ithe above procedure, we ap- 
ply it to the sillium modeD, one of the prototype mod- 
els to study network- forming materials. In this model, 
tuned for amorphous silicon, an explicit list of covalent 
bonds between pairs of Si atoms is kept, with the prop- 
erty that each Si atom is bonded to four neighboring 
atoms. The energy is described by the Keating poten- 
tial, which contains a quadratic penalty for bond-length 
deviations from the crystalline distance of 2.35 A, and 
a quadratic penalty for bond-angle deviations from the 
tetrahedral angle 8o = arccos(— 1/3). The list of bonds 
determines the atomic positions uniquely, since in this 
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model the energy is minimized at all times. As a result, 
the phase space of this model is limited to a finite number 
of 3A^-dimensional points. The evolution of the network 
consists of a large number of random bond transpositions, 
each accepted with the Metropolis probability: 



P = min 



1, exp 



(5) 



where T is the temperature and and Ef are the total 
(minimized) energies of the system before and after the 
proposed bond transposition. A typical sillium network 
is shown in Fig. 3- Our networks are generated as de- 
scribed in Refs. oH, in which a number of algorithmic 
improvements have been proposed as compared to the 
original algorithm of Wooten, Winer and Weaire. 

To verify the validity of the information theory ap- 
proach we first turn to a two-dimensional version of the 
sillium model. Contrary to the three-dimensional sillium 
model, the dynamics here does not lead to glassy states. 
This allows us to determine the entropy in equilibrium us- 
ing both the standard approach — via Eq. (|l|) — as well 
as information theory. In this two-dimensional model, 
atoms are three-fold coordinated and the ideal bond an- 
gle is 120 degrees; a typical configuration is shown in 
Fig. I 

We simulate a two-dimensional sillium network con- 
taining 1008 atoms for a number of (physically interest- 
ing) temperatures ranging from ksT = 0.15 eV to ksT = 
0.60 eV using periodic boundary conditions in both di- 
rections. At each temperature we bring the system to 
equilibrium with 100 attempted bond-transpositions per 
atom. Next, 20 snapshots of the network are stored; 
each separated by five attempted bond transpositions per 
atom. From these snapshots the average energy and the 
Shannon entropy are obtained. In determining the Shan- 
non entropy we draw clusters from each snapshot simul- 
taneously. After the simulation has covered the entire 
temperature range, the entropy is also obtained from the 
average energy measurements and Eq. (^. 

In Fig. ^ we show the configurational entropy as a func- 
tion of temperature. The solid curve shows the entropy 
per atom as obtained using the standard thermodynamic 
approach; crosses show the entropy per atom as obtained 
using information theory. There appears to be a phase 
transition near ksT = 0.37 eV shown by the sudden 
jump in entropy. The remarkable feature of the displayed 
results is the overall good agreement between the two 
methods, even near the phase transition. 

Next, we use information theory to determine the con- 
figurational entropy of a well-relaxed three-dimensional 
sillium network containing 20,000 atoms. This network 
was generated using the improved Wooten, Winer and 
Weaire algorithmuQ. The Keating energy of this net- 
work is 0.286 eV per atom; the standard deviation in the 
mean bond-angle is 9.63 degrees. Structural and elec- 
tronic properties of this network are in excellent agree- 
ment with experimental properties of amorphous siliconB. 
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FIG. 4: Diagram showing the configurational entropy per 
atom in units of fcs as a function of temperature in eV for 
a two-dimensional sillium network consisting of 1008 atoms. 
The solid curve shows the entropy per atom as obtained us- 
ing the standard thermodynamic approach; crosses show the 
entropy per atom obtained using information theory. 



Applying our method to this model, we obtain a config- 
urational entropy of l.02kB per atom. 

The (classical) vibrational entropy is obtained from the 
eigenvalues of the dynamical matrix and for this network 
was found to be 3. Oiks per atom. For the crystalline 
phase the configurational entropy is zero and only the vi- 
brational entropy contributes to the entropy, in this case 
3.18fcB per atom. The difference in entropy between the 
crystalline and amorphous phases of silicon is therefore 
estimated to be 0.91A:b per atom. 

The energy difference between the crystalline and 
amorphous phases of silicon was recently determined by 
BiswasQ using tight-binding, who found an energy differ- 
ence of 0.18 eV per atom. 

With the above estimates for the differences in energy 
and entropy we can calculate the transition temperature 
Tc, defined as the temperature where the free energy dif- 
ference between the crystalline phase and the amorphous 
phase changes sign. It is given by: 



Tr. = 



AE 
'AS' 



(6) 



where AE and AS are the energy difference and the en- 
tropy difference between the crystalline and the amor- 
phous phase, respectively. Substitution of our estimates 
for AE and AS* into Eq. (|) yields « 2300K. This 
temperature compares well with the value of Tc « 2500K 
as inferred from calorimetric experimentsu, and is well 
above the melting point of silicon. It confirms that the 
amorphous phase is not thermodynamically stable at any 
temperature. 

As a final application of our method we determine the 
configurational entropy of vitreous silica. The structure 
of this material is formed by covalent bonds between sili- 
con atoms and oxygen atoms. Barring rare defects, each 
silicon atom is bonded to four oxygen atoms and each 
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oxygen atom is bonded to two silicon atoms. We gen- 
erate a silica network containing 3000 atoms (with pe- 
riodic boundary conditions) in a spirit similar to amor- 
phous silicon, but we replace the Keating potential by 
the one of Tu and TersoflHlMj. The resulting network-is 
then quenched with the more realistic BKS potentialtil, 
using parameters as described in Ref. |l2[ The BKS en- 
ergy per silicon atom of the quenched network is 0.13 eV 
above that of the /3-cristobalite structure. For compar- 
ison, samples prepared by molecular dynamics typically 
yield a much larger energy difference of 0.30 eV. This 
clearly demonstrates that our network is well relaxed. 
After quenching with the BKS potential, the atoms re- 
main perfectly coordinated: the silicon-oxygen bonds are 
easily reconstructed from a distance criterion and the 
silicon-oxygen radial distribution function. For our net- 
work, the average 0-Si-O bond-angle is 109.43 degrees 
with a standard deviation of 4.40 degrees; the average 
Si-O-Si bond-angle was found to be 150.85 degrees with 
a standard deviation of 12.09 degrees. 

To estimate the configurational entropy per silicon 
atom in well-relaxed vitreous silica, we used the above 



network and replaced each oxygen atom plus its two 
bonds by a single silicon-silicon bond. Next, we ap- 
plied the information-theoretic method described earlier 
on this network and found for the configurational en- 
tropy 0.97^5 per silicon atom. Given the limited size of 
this sample we estimate that the actual entropy might be 
up to 10% larger. Ignoring the difference in vibrational 
entropy, as well as the entropy contribution of possible 
multiple oxygen positions, we obtain for the transition 
temperature of vitreous silica ~ 1600K. 

In summary, we have developed a computationally effi- 
cient method to determine the configurational entropy of 
network forming materials. For well-relaxed samples of 
amorphous silicon and vitreous silica, we find for the en- 
tropy per silicon atom 1.02/cb and 0.97A:b, respectively. 
In future research, we hope to extend the applicability 
of this method to other disordered materials such as col- 
loidal systems and metallic glasses. 

We thank Henk van Beijeren and Normand Mousseau 
for stimulating discussion, and Normand Mousseau also 
for the computation of the BKS energies. 
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